import numpy as np
import matplotlib.pyplot as plt
# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False    # 用来正常显示负号

# 设置参数
a = 1.0  # x方向长度
b = 1.0  # y方向长度
n_terms = 10  # 级数截断项数

# 创建网格
x = np.linspace(0, a, 200)
y = np.linspace(0, b, 200)
X, Y = np.meshgrid(x, y)

# 初始化解
U = X * (a - X)  # 特解部分 -x(x-a) = x(a-x)

# 计算级数部分
for k in range(1, n_terms + 1):
    n = 2 * k - 1  # (2k-1)
    denominator = n**3 * np.pi**3 * np.cosh(n * np.pi * b / (2 * a))
    coeff = 8 * a**2 / denominator
    
    # 计算级数项
    sinh_arg = n * np.pi * (Y - b/2) / a
    sin_arg = n * np.pi * X / a
    term = coeff * np.cosh(sinh_arg) * np.sin(sin_arg)
    
    U -= term  # 从解中减去级数项

# 创建二维图像
plt.figure(figsize=(10, 8))

# 绘制伪彩色图
contour = plt.contourf(X, Y, U, cmap='viridis', levels=50)
plt.colorbar(contour, label='u(x,y)')

# 添加等高线
C = plt.contour(X, Y, U, colors='black', linewidths=0.5, levels=8)
plt.clabel(C, inline=True, fontsize=8, fmt='%1.3f')

# 添加标题和标签
plt.title('泊松方程解 $\Delta u = -2$ (边界条件为0)', fontsize=14)
plt.xlabel('x', fontsize=12)
plt.ylabel('y', fontsize=12)

# 添加网格
plt.grid(True, linestyle='--', alpha=0.3)

# 标记中心点
plt.plot(0.5, 0.5, 'ro', markersize=5)
plt.annotate('最大值 (0.073)', xy=(0.5, 0.5), xytext=(0.6, 0.6),
             arrowprops=dict(facecolor='red', shrink=0.05, width=1))

# 设置坐标轴范围
plt.xlim(0, a)
plt.ylim(0, b)

# 以300dpi保存图像
plt.savefig('例题8-2.png', dpi=300, bbox_inches='tight')
plt.close()

print("二维平面图已保存为 '例题8-4.png' (300dpi)")